import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint

plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
plt.rc('font', size=10)
plt.rc('font', family='SimHei')

m_A = 4866
m_B = 2433
f = 1760
delta_m = 1091.099
omiga = 1.9802
k_3 = 528.5018  # 兴波
k_5 = 80000  # 弹簧劲度系数
k_6 = 10000  # 阻尼系数
ro = 1025
g = 9.8
s = np.pi
l_0 = 0.5  # 弹簧原长
I_1 = 8289.43  # 浮子转动惯量
h_2 = 0.5  # 振子的高度
delta_I = 7001.914
k_7 = 1655.909  # 兴波
k_8 = 250000  # 弹簧劲度系数
k_9 = 1000  # 阻尼系数
k_10 = 8890.7  # 静水恢复力矩
L = 2140

# 在SA函数外部初始化记录容器
history = {
    'temperature': [],
    'best_P': [],
    'current_P': [],
    'k6': [],
    'k9': [],
}


def plot_SA_history():
    plt.figure(figsize=(15, 8))

    # 子图1：损失函数变化
    plt.subplot(2, 2, 1)
    plt.plot(history['best_P'], 'r-', label='Best P')
    plt.plot(history['current_P'], 'b--', alpha=0.5, label='Current P')
    plt.xlabel('Iteration')
    plt.ylabel('P')
    plt.title('P变化曲线')
    plt.legend()
    plt.grid(True, which="both", ls="--")

    # 子图2：温度衰减曲线
    plt.subplot(2, 2, 2)
    plt.plot(history['temperature'], 'g-', label='temperture')
    plt.xlabel('Iteration')
    plt.ylabel('Temperature')
    plt.title('温度下降曲线')
    plt.grid(True, ls="--")

    # 子图3：k6参数演化
    plt.subplot(2, 2, 3)
    plt.plot(history['k6'])
    plt.xlabel('Iteration')
    plt.ylabel('k6参数值')
    plt.title('k6参数变化曲线')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, ls="--")

    # 子图3：k9参数演化
    plt.subplot(2, 2, 4)
    plt.plot(history['k9'])
    plt.xlabel('Iteration')
    plt.ylabel('k9参数值')
    plt.title('k9参数变化曲线')
    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
    plt.grid(True, ls="--")

    plt.tight_layout()
    plt.show()


def P(k_6, k_9):
    # 定义一个方程组（微分方程组）
    # 计算振子的转动惯量
    def I_2(x_1, x_2):
        delta_x = m_B * g / k_5
        l_1 = x_2 - x_1 + l_0 - delta_x
        l_2 = l_1 + h_2
        return 1 / 3 * m_B / h_2 * (l_2 ** 3 - l_1 ** 3)

    # 定义一个方程组（微分方程组）
    def pfun(y, t):
        y0, y1, y2, y3, y4, y5, y6, y7 = y
        return np.array([
            1 / (m_A + delta_m) * (
                        -k_3 * y0 - ro * g * s * y1 + f * np.cos(omiga * t) + k_5 * (y3 - y1) + k_6 * (y2 - y0)),
            y0,
            1 / m_B * (-k_5 * (y3 - y1) - k_6 * (y2 - y0)),
            y2,
            1 / (I_1 + delta_I) * (-k_7 * y4 - k_10 * y5 + L * np.cos(omiga * t) + k_8 * (y7 - y5) + k_9 * (y6 - y4)),
            y4,
            1 / I_2(y1, y3) * (-k_8 * (y7 - y5) - k_9 * (y6 - y4)),
            y6,
        ])

    t = np.arange(0, 1000, 0.2)  # 创建自变量序列
    soli = odeint(pfun, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], t)  # 求数值解

    re = 0
    for i in range(200, len(soli[:, 2])):
        re += (k_6 * (soli[:, 0][i] - soli[:, 2][i]) ** 2 + k_9 * (soli[:, 6][i] - soli[:, 4][i]) ** 2) * 0.2

    return re / 960


# 自适应模拟退火算法
def adaptive_SA(iter, t0, tf, alpha_base):
    # 初始化自适应参数
    t = t0
    k6c = 20000
    k9c = 20000
    Pc = P(k6c, k9c)
    k6b = k6c
    k9b = k9c
    Pb = Pc

    # 自适应参数
    step_size = 50000  # 初始扰动范围
    accept_rate = 0.0  # 接受率
    accept_count = 0  # 接受次数
    no_improve_count = 0  # 未改进计数
    alpha = alpha_base

    for i in range(iter):
        # 自适应调整扰动范围
        if i > 50:  # 等待初始阶段结束
            if accept_rate > 0.5:
                step_size *= 1.1  # 接受率高，增大搜索范围
            else:
                step_size *= 0.9  # 接受率低，减小搜索范围
            step_size = np.clip(step_size, 1000, 100000)  # 限制扰动范围

        # 生成新解 - 使用当前扰动范围
        k6n = k6c + np.random.randint(-step_size, step_size)
        k6n = np.clip(k6n, 0, 100000)
        k9n = k9c + np.random.randint(-step_size, step_size)
        k9n = np.clip(k9n, 0, 100000)

        Pn = P(k6n, k9n)

        # 自适应接受概率
        delta = Pn - Pc
        accept_prob = np.exp(delta / t) if delta < 0 else 1.0

        # 决定是否接受新解
        if Pn > Pc or np.random.rand() < accept_prob:
            k6c = k6n
            k9c = k9n
            Pc = Pn
            accept_count += 1

            # 更新最优解
            if Pc > Pb:
                k6b = k6c
                k9b = k9c
                Pb = Pc
                no_improve_count = 0  # 重置未改进计数
            else:
                no_improve_count += 1
        else:
            no_improve_count += 1

        # 计算接受率 (基于最近50次迭代)
        if i >= 50:
            accept_rate = accept_count / 50
            accept_count = 0  # 重置计数器

        # 自适应降温策略
        if no_improve_count > 20:  # 长期未改进
            alpha = min(alpha_base * 1.1, 0.99)  # 加速降温
        else:
            alpha = alpha_base  # 恢复基础降温速率

        t = t * alpha

        # 提前终止条件
        if t < tf or (i > 100 and no_improve_count > 50):
            break

        print(
            f"第{i}轮, 温度:{t:.2f}, 最优P:{Pb:.2f}, 当前P:{Pc:.2f}, 扰动范围:{step_size:.0f}, 接受率:{accept_rate:.2f}")

        # 记录历史数据
        history['temperature'].append(t)
        history['best_P'].append(Pb)
        history['current_P'].append(Pc)
        history['k6'].append(k6c)
        history['k9'].append(k9c)

    return k6b, k9b, Pb


re = adaptive_SA(800, 2000, 0.001, 0.9)
print(re)
plot_SA_history()

for key in history.keys():
    history[key].clear()
print(history['temperature'])
